Librerías:

library(tidyverse)
library(data.table)
library(readxl)
library(GGally)
library(corrplot)
library(leaflet)
library(raster)
library(mltools)
library(highcharter)

Lectura de datos:

raw_encuesta <- fread("databases/calidad_vida_ok.csv", encoding = "UTF-8") %>% 
  as_tibble()

# Unificación de las estadísticas de los barrios "ALTAVISTA" y "CABECERA ALTAVISTA"

raw_encuesta <- raw_encuesta %>% mutate(
  encuesta_calidad.barrio = case_when(
    encuesta_calidad.barrio == "CABECERA ALTAVISTA" ~ "ALTAVISTA",
    encuesta_calidad.barrio == "CIUDADELA NUEVO OCCIDENTE" ~ "CABECERA URBANA CORREGIMIENTO SAN CRISTÓBAL",
    encuesta_calidad.barrio == "AREA DE EXPANCION SAN CRISTOBAL" ~ "ÁREA DE EXPANSIÓN SAN CRISTÓBAL",
    encuesta_calidad.barrio == "CABECERA SAN CRISTÓBAL" ~ "CABECERA URBANA CORREGIMIENTO SAN CRISTÓBAL",
    encuesta_calidad.barrio == "SAN CRISTOBAL" ~ "CABECERA URBANA CORREGIMIENTO SAN CRISTÓBAL",
    encuesta_calidad.barrio == "PROGRESO  Nº 2" ~ "PROGRESO",
    TRUE ~ encuesta_calidad.barrio
  )
)

diccionario <- read_excel("databases/diccionario_ecv_2016.xlsx")

Número de encuestados:

raw_encuesta %>% 
  count(`encuesta_calidad.barrio`, sort = T) %>% 
  arrange(n)
## # A tibble: 296 x 2
##    encuesta_calidad.barrio     n
##    <chr>                   <int>
##  1 EL CARMELO                  4
##  2 GUAYAQUIL                   4
##  3 LA ILUSIÓN                  6
##  4 SAN JOSÉ                   10
##  5 BOQUERÓN                   15
##  6 EL PATIO                   19
##  7 BARRIO COLOMBIA            21
##  8 DESCONOCIDO                21
##  9 YARUMALITO                 22
## 10 EL UVITO                   25
## # ... with 286 more rows

Creación de la base de datos:

encuesta_ingresos<-raw_encuesta%>%
  mutate(proporcion_ingreso_extra=case_when(encuesta_calidad.p_71 == 1 ~1,
                                 TRUE~ 0))%>%
  mutate(proporcion_pension=case_when(encuesta_calidad.p_71 != 0~1,
                                 TRUE~ 0))%>%
  mutate(subsidio=case_when(encuesta_calidad.p_93 != 0~1,
                                 TRUE~ 0))%>%
  mutate(salario=case_when(encuesta_calidad.p_121 != 0~1,
                                 TRUE~ 0))%>%
  mutate(arriendo=case_when(encuesta_calidad.p_227 != 0~1,
                                 TRUE~ 0))%>%
  mutate(total_gastos=case_when(encuesta_calidad.p_232 != 0~1,
                                 TRUE~ 0))%>%
  mutate(gasto_servicios=case_when(encuesta_calidad.p_240 != 0~1,
                                 TRUE~ 0))%>%
  mutate(gasto_hogar=case_when(encuesta_calidad.p_254 != 0~1,
                                 TRUE~ 0))%>%
  mutate(encuesta_calidad.barrio = str_replace(encuesta_calidad.barrio, "ANDALUCIA", "ANDALUCÍA") %>% 
  str_replace("Nº 2", "NO.2") %>% 
  str_replace("Nº 1", "NO.1") %>% 
  str_replace("Nº 3", "NO.3") %>%
  str_replace("AREA EXPANSION", "ÁREA DE EXPANSIÓN") %>%
  str_replace("EXPANCION", "EXPANSIÓN") %>% 
  str_replace("AREA", "ÁREA") %>% 
  str_replace("BOMBONA", "BOMBONÁ") %>% 
  str_replace("LA ASOMADERA", "ASOMADERA") %>%
  str_replace("BELALCAZAR", "BELALCÁZAR") %>% 
  str_replace("CALAZANS", "CALASANZ") %>% 
  str_replace("COLON", "COLÓN") %>% 
  str_replace("MIRA FLORES", "MIRAFLORES") %>% 
  str_replace("BARRIO FACULTAD DE MINAS", "FACULTAD DE MINAS") %>% 
  str_replace("CABECERA SAN ANT DE PR.", "SAN ANTONIO DE PRADO") %>% 
  str_replace("CARLOS E RESTREPO", "CARLOS E. RESTREPO") %>% 
  str_replace("URQUITA", "URQUITÁ") %>% 
  str_replace("LOS CERROS EL VERJEL", "LOS CERROS EL VERGEL") %>% 
  str_replace("CAYCEDO", "CAICEDO") %>% 
  str_replace("VALDES", "VALDÉS") %>% 
  str_replace("CERRO EL VOLADOR", "B. CERRO EL VOLADOR") %>% 
  str_replace("MOSCU", "MOSCÚ") %>% 
  str_replace("JOSELA", "JOSÉ LA") %>%
  str_replace("JOSE", "JOSÉ") %>% 
  str_replace("EL YOLOMBO", "YOLOMBO") %>% 
  str_replace("PIEDRAS BLANCAS", "PIEDRAS BLANCAS - MATASANO") %>% 
  str_replace("BASILIA", "BRASILIA") %>% 
  str_replace("VILLA TINA", "VILLATINA") %>% 
  str_replace("LILIAM", "LILLIAM") %>% 
  str_replace("BOLIVAR", "BOLÍVAR") %>% 
  str_replace("CORREGIMIENTO PALMITAS", "PALMITAS SECTOR CENTRAL") %>% 
  str_replace("INES", "INÉS") %>% 
  str_replace("FE", "FÉ") %>% 
  str_replace("LUCIA", "LUCÍA") %>% 
  str_replace("SABIO", "SAVIO") %>% 
  str_replace("BERMEJAL- LOS ÁLAMOS", "BERMEJAL-LOS ÁLAMOS") %>% 
  str_replace("BOLÍVARIANA", "BOLIVARIANA") %>% 
  str_replace("EL NOGAL - LOS ALMENDROS", "EL NOGAL-LOS ALMENDROS") %>% 
  str_replace("JUAN XXIII - LA QUIEBRA", "JUAN XXIII LA QUIEBRA") %>% 
  str_replace("PROGRESO  Nº 2", "EL PROGRESO") %>% 
  str_replace("MARIA", "MARÍA") %>% 
  str_replace("PLAYÓN", "PLAYON") %>% 
  str_replace("EL SOCORRO / LA GABRIELA", "EL SOCORRO") %>% 
  str_replace("FÉRRINI", "FERRINI") %>% 
  str_replace("LA CANDE LARIA", "LA CANDELARIA") %>%
  str_replace("EL PLAYON", "PLAYÓN") %>%
  str_replace("IGUANA", "IGUANÁ") %>%
  str_replace("MARÍA CANO - CARAMBOLAS", "MARÍA CANO-CARAMBOLAS") %>%
  str_replace("DE ABURRA", "DEL ABURRÁ") %>%
  str_replace("ALTAVISTA CENTRAL", "ALTAVISTA SECTOR CENTRAL") %>%
  str_replace("SECTOR CENTRAL", "CENTRO ADMINISTRATIVO") %>%
  str_replace("ALTAVISTA CENTRO ADMINISTRATIVO", "ALTAVISTA SECTOR CENTRAL") %>%
  str_replace("SANTA ELENA CENTRO ADMINISTRATIVO", "SANTA ELENA SECTOR CENTRAL") %>%
  str_replace("PALMITAS CENTRO ADMINISTRATIVO", "PALMITAS SECTOR CENTRAL") %>%  
  str_replace("PROGRESO", "EL PROGRESO")
  )



db_ingresos<-encuesta_ingresos%>%
  group_by(encuesta_calidad.barrio, encuesta_calidad.comuna)%>%
  summarize(n=n(),
            proporcion_ingreso_extra=sum(proporcion_ingreso_extra==1,na.rm=TRUE)/sum(!is.na(proporcion_ingreso_extra),na.rm=TRUE),
            proporcion_pension=sum(proporcion_pension==1,na.rm = TRUE)/sum(!is.na(proporcion_pension),na.rm=TRUE),
            subsidio=sum(`encuesta_calidad.p_93`, na.rm=TRUE)/sum(subsidio==1,na.rm = TRUE),
            salario=sum(`encuesta_calidad.p_121`, na.rm=TRUE)/sum(salario==1,na.rm = TRUE),
            arriendo=sum(`encuesta_calidad.p_227`,na.rm=TRUE)/sum(arriendo==1,na.rm = TRUE),
            total_gastos=sum(`encuesta_calidad.p_232`,na.rm=TRUE)/sum(total_gastos==1,na.rm = TRUE),
            gasto_servicios=sum(`encuesta_calidad.p_240`,na.rm=TRUE)/sum(gasto_servicios==1,na.rm = TRUE)
            )%>%ungroup()%>%
  mutate(subsidio= replace_na(subsidio, 0))%>%
  mutate(salario= replace_na(salario, 0))%>%
  mutate(arriendo= replace_na(arriendo, 0))%>%
  mutate(total_gastos= replace_na(total_gastos, 0))%>%
  mutate(gasto_servicios= replace_na(gasto_servicios, 0))

# Eliminamos el barrio Desconocido
db_ingresos <- db_ingresos %>% 
  filter(encuesta_calidad.barrio != "DESCONOCIDO")

db_ingresos[complete.cases(db_ingresos),]
## # A tibble: 308 x 10
##    encuesta_calida~ encuesta_calida~     n proporcion_ingr~
##    <chr>            <chr>            <int>            <dbl>
##  1 AGUAS FRÍAS      ALTAVISTA          126          0.00794
##  2 ALDEA PABLO VI   POPULAR            426          0.00235
##  3 ALEJANDRÍA       EL POBLADO         709          0.00987
##  4 ALEJANDRO ECHAV~ BUENOS AIRES      1094          0.00640
##  5 ALFONSO LÓPEZ    CASTILLA          1994          0.00752
##  6 ALTAMIRA         ROBLEDO            755          0.00927
##  7 ALTAVISTA        ALTAVISTA          168          0.0119 
##  8 ALTAVISTA        BELEN             1988          0.00553
##  9 ALTAVISTA SECTO~ ALTAVISTA          507          0.00592
## 10 ALTOS DEL POBLA~ EL POBLADO         375          0.00267
## # ... with 298 more rows, and 6 more variables: proporcion_pension <dbl>,
## #   subsidio <dbl>, salario <dbl>, arriendo <dbl>, total_gastos <dbl>,
## #   gasto_servicios <dbl>

Escritura.

write_excel_csv2(db_ingresos, "databases/db_ingresos.csv")
write_excel_csv2(db_ingresos[complete.cases(db_ingresos),], "databases/db_ingresos_completecases.csv")

Esta dimensión llamada “ingresos” aborda algunas preguntas referentes a los ingresos y gastos de las personas de los diferentes barrios de Medellín.

DESCRIPCIÓN DE LAS VARIABLES:

Análisis descriptivo

summary(db_ingresos)
##  encuesta_calidad.barrio encuesta_calidad.comuna       n         
##  Length:308              Length:308              Min.   :   4.0  
##  Class :character        Class :character        1st Qu.: 348.8  
##  Mode  :character        Mode  :character        Median : 862.5  
##                                                  Mean   :1073.2  
##                                                  3rd Qu.:1439.0  
##                                                  Max.   :8987.0  
##  proporcion_ingreso_extra proporcion_pension    subsidio      
##  Min.   :0.000000         Min.   :0.1000     Min.   :      0  
##  1st Qu.:0.003061         1st Qu.:0.4049     1st Qu.: 244619  
##  Median :0.005511         Median :0.4426     Median : 370100  
##  Mean   :0.006463         Mean   :0.4386     Mean   : 520659  
##  3rd Qu.:0.008005         3rd Qu.:0.4696     3rd Qu.: 587997  
##  Max.   :0.066667         Max.   :0.7500     Max.   :8241967  
##     salario           arriendo        total_gastos     gasto_servicios 
##  Min.   :      0   Min.   :      0   Min.   : 425763   Min.   : 50500  
##  1st Qu.: 305632   1st Qu.: 283998   1st Qu.: 833227   1st Qu.:150364  
##  Median : 437404   Median : 409747   Median :1087169   Median :206197  
##  Mean   : 703446   Mean   : 535185   Mean   :1393247   Mean   :225134  
##  3rd Qu.: 742865   3rd Qu.: 654791   3rd Qu.:1580345   3rd Qu.:275162  
##  Max.   :8734250   Max.   :1901874   Max.   :4976319   Max.   :588615
min_gasto_servicios<-apply(db_ingresos[10],2,min, na.rm=T)
db_ingresos[which(db_ingresos$gasto_servicios==min_gasto_servicios),1]
## # A tibble: 1 x 1
##   encuesta_calidad.barrio
##   <chr>                  
## 1 SAN JOSÉ
max_gasto_servicios<-apply(db_ingresos[10],2,max, na.rm=T)
db_ingresos[which(db_ingresos$gasto_servicios==max_gasto_servicios),1]
## # A tibble: 1 x 1
##   encuesta_calidad.barrio
##   <chr>                  
## 1 ALTOS DEL POBLADO

Gráfico (1).

grafico<-ggpairs(db_ingresos,columns = 4:10, mapping=aes(alpha = 0.4))
ggsave(filename="images/correlacion.png", plot =grafico,  width = 20, height = 20)
Correlación entre las variables

Correlación entre las variables

Matriz de correlaciones:

data<-as.matrix(db_ingresos[,3:9])
mcor<-cor(data)
mcor
##                                    n proporcion_ingreso_extra
## n                         1.00000000              -0.06518933
## proporcion_ingreso_extra -0.06518933               1.00000000
## proporcion_pension       -0.17011586              -0.03081172
## subsidio                 -0.01152908              -0.02913350
## salario                  -0.02949656              -0.01847511
## arriendo                 -0.10414827               0.01745026
## total_gastos             -0.10706213               0.01250525
##                          proporcion_pension    subsidio     salario
## n                               -0.17011586 -0.01152908 -0.02949656
## proporcion_ingreso_extra        -0.03081172 -0.02913350 -0.01847511
## proporcion_pension               1.00000000  0.03797035  0.08478863
## subsidio                         0.03797035  1.00000000  0.28460736
## salario                          0.08478863  0.28460736  1.00000000
## arriendo                         0.05240449  0.40708611  0.65315248
## total_gastos                     0.02899374  0.40252794  0.63576569
##                             arriendo total_gastos
## n                        -0.10414827  -0.10706213
## proporcion_ingreso_extra  0.01745026   0.01250525
## proporcion_pension        0.05240449   0.02899374
## subsidio                  0.40708611   0.40252794
## salario                   0.65315248   0.63576569
## arriendo                  1.00000000   0.95678984
## total_gastos              0.95678984   1.00000000

Gráfico (2) Correlaciones:

corrplot(mcor,type = "upper",order = "hclust",tl.col = "black",tl.srt = 45)

Escalemos.

db_ingresos_escalado <- db_ingresos%>% 
  mutate_at(vars(-n, -encuesta_calidad.barrio, -encuesta_calidad.comuna), scale)

Selección de las columnas para el agrupamiento.

datos_escalados_agrupamiento <- db_ingresos_escalado %>% 
  dplyr::select(-n, -encuesta_calidad.barrio, -encuesta_calidad.comuna)

Apliquemos análisis de componente principales (PCA).

pca <- prcomp(datos_escalados_agrupamiento)
summary(pca)
## Importance of components:
##                           PC1    PC2    PC3    PC4     PC5    PC6     PC7
## Standard deviation     1.8902 1.0168 0.9845 0.8853 0.68540 0.3579 0.20592
## Proportion of Variance 0.5104 0.1477 0.1385 0.1120 0.06711 0.0183 0.00606
## Cumulative Proportion  0.5104 0.6581 0.7966 0.9085 0.97564 0.9939 1.00000

Varianza explicada por cada componente

porcentaje_exp <- 100 * (pca$sdev ^ 2) / sum(pca$sdev ^ 2)
ggplot(data = data.frame(componente = factor(1:ncol(datos_escalados_agrupamiento)), porcentaje_exp)) +
  geom_col(mapping = aes(x = componente, y = porcentaje_exp))

suma_acumulada_porcentaje <- cumsum(porcentaje_exp)
ggplot(data = data.frame(componente = factor(1:ncol(datos_escalados_agrupamiento)), suma_acumulada_porcentaje),
        aes(x = componente, y = suma_acumulada_porcentaje, group = 1)) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 90, color = "red") +
  geom_vline(xintercept = 7)

suma_acumulada_porcentaje
## [1]  51.03878  65.80981  79.65645  90.85302  97.56407  99.39425 100.00000

Proyección respecto a las componentes principales

datos_proyectados <- t(t(pca$rotation) %*% t(datos_escalados_agrupamiento))
datos_agrupamiento <- datos_proyectados[, 1:4] %>% as_tibble()

Agrupamiento mediante K-means

Se utilizará la metodología de Elbow para identificar el k óptimo que ofreza una menor distancia entre clusters como medida de error Wss.

datos_escalados_agrupamiento<- db_ingresos_escalado %>% 
  dplyr::select(-n, -encuesta_calidad.barrio, -encuesta_calidad.comuna)
k_maximo <- 15
wss <- sapply(1:k_maximo, 
              function(k) {
                kmeans(datos_escalados_agrupamiento, k, nstart = 10, iter.max = 15 )$tot.withinss
              })
plot(1:k_maximo,wss,
     type="b", pch = 19, frame = FALSE, 
     xlab="Número de grupos [K]",
     ylab="Total de suma de cuadrados intra-grupos")

Según nuestro algoritmo de K-means, nuestro K óptimo es de 9 grupos, pero se considera más pertinente ultilizar un número menor para favorecer la interpretabilidad.

k_optimo_kmeans<-9
k_seleccionado<-5

Agrupamiento con K óptimo

Fijaremos una semilla de los número aleatorias para garantizar la reproducibilidad de las conclusiones del agrupamiento.

set.seed(3)
agrupamiento_kmeans <- kmeans(datos_agrupamiento, centers = k_seleccionado, nstart = 10, iter.max = 15 )
db_ingresos_kmeans <- bind_cols(db_ingresos, cluster = agrupamiento_kmeans$cluster)

Mapa geográfico mediante el agrupamiento K-means

political <- shapefile("Barrio_Vereda/Barrio_Vereda.shp")
Encoding(political@data$NOMBRE) <- "UTF-8"
political$NOMBRE <- political$NOMBRE %>% toupper() %>% str_replace("DE  MESA", "DE MESA")

grupos_barrios <- data.frame(barrio_nombre = db_ingresos$encuesta_calidad.barrio, grupo = agrupamiento_kmeans$cluster)
nombres_mapa <- data.frame(nombre_barrio = political$NOMBRE)

vector_nombres = c()
vector_grupos = c()

for(nombre_mapa in nombres_mapa$nombre_barrio) {
  grupo <- grupos_barrios[grupos_barrios$barrio_nombre == nombre_mapa, 2][1]
  vector_nombres <- c(vector_nombres, nombre_mapa)
  vector_grupos <- c(vector_grupos, grupo)
}


factpal <- colorFactor(rainbow(k_seleccionado), vector_grupos)

mapa_grupos <-tibble(
  nombre_barrio = vector_nombres, 
  grupo = vector_grupos,
  color = factpal(grupo),
  descripcion = case_when(
    grupo == 1 ~ "Grupo 1",
    grupo == 2 ~ "Grupo 2",
    grupo == 3 ~ "Grupo 3",
    grupo == 4 ~ "Grupo 4",
    grupo == 5 ~ "Grupo 5",
    TRUE ~ "Sin muestra representativa")
)

colores <- mapa_grupos %>%
  dplyr::select(-nombre_barrio) %>% 
  distinct(color, .keep_all = TRUE) %>% 
  arrange(desc(grupo))


leaflet(data = political) %>% 
  addTiles() %>% 
  addProviderTiles(providers$OpenStreetMap) %>% 
  addPolygons(fill = TRUE, stroke = TRUE, weight = 2, color = mapa_grupos$color, 
              label = as.character(political$NOMBRE),
              popup = as.character(political$NOMBRE)) %>% 
  addLegend("bottomright", colors = colores$color, labels = as.character(colores$descripcion))

Análisis descriptivo de los grupos obtenido con K-means

grupo_1<-grupos_barrios[grupos_barrios[,2]==1,]
grupo_2<-grupos_barrios[grupos_barrios[,2]==2,]
grupo_3<-grupos_barrios[grupos_barrios[,2]==3,]
grupo_4<-grupos_barrios[grupos_barrios[,2]==4,]
grupo_5<-grupos_barrios[grupos_barrios[,2]==5,]

Tamaño del grupo:

nrow(grupo_1)
## [1] 9
nrow(grupo_2)
## [1] 29
nrow(grupo_3)
## [1] 100
nrow(grupo_4)
## [1] 74
nrow(grupo_5)
## [1] 96

Para este análisis se utilizará la siguiente convención de medidas descriptivas: media (desviación estándar).**

media_sd <- function(variable) {
  paste0(round(mean(variable), 3), " (", round(sd(variable), 3), ")")
}

db_ingresos_kmeans %>%
  dplyr::select(-encuesta_calidad.barrio, -encuesta_calidad.comuna, -n) %>% 
  group_by(cluster) %>% 
  summarise_all(~media_sd(.))
## # A tibble: 5 x 8
##   cluster proporcion_ingr~ proporcion_pens~ subsidio salario arriendo
##     <int> <chr>            <chr>            <chr>    <chr>   <chr>   
## 1       1 0.033 (0.013)    0.419 (0.061)    300988.~ 180944~ 304928.~
## 2       2 0.008 (0.006)    0.438 (0.042)    1511526~ 248279~ 1373217~
## 3       3 0.006 (0.004)    0.473 (0.042)    309749.~ 337484~ 318650.~
## 4       4 0.006 (0.003)    0.459 (0.031)    657956.~ 100480~ 762689.~
## 5       5 0.005 (0.003)    0.389 (0.037)    355791.~ 363832~ 353805.~
## # ... with 2 more variables: total_gastos <chr>, gasto_servicios <chr>

Agrupamiento jeráquico:

cantidad_grupos_jerarquico <- 5

distancia <- dist(datos_agrupamiento)
agrupamiento <- hclust(distancia, method = "complete")
arbol_corte <- cutree(agrupamiento, k = cantidad_grupos_jerarquico)

Mapa geográfico usando el agrupamiento Jerárquico

political <- shapefile("Barrio_Vereda/Barrio_Vereda.shp")
Encoding(political@data$NOMBRE) <- "UTF-8"
political$NOMBRE <- political$NOMBRE %>% toupper() %>% str_replace("DE  MESA", "DE MESA")

grupos_barrios <- data.frame(barrio_nombre = db_ingresos$encuesta_calidad.barrio, grupo =arbol_corte)
nombres_mapa <- data.frame(nombre_barrio = political$NOMBRE)

vector_nombres = c()
vector_grupos = c()

for(nombre_mapa in nombres_mapa$nombre_barrio) {
  grupo <- grupos_barrios[grupos_barrios$barrio_nombre == nombre_mapa, 2][1]
  vector_nombres <- c(vector_nombres, nombre_mapa)
  vector_grupos <- c(vector_grupos, grupo)
}


factpal <- colorFactor(rainbow(k_seleccionado), vector_grupos)

mapa_grupos <-tibble(
  nombre_barrio = vector_nombres, 
  grupo = vector_grupos,
  color = factpal(grupo),
  descripcion = case_when(
    grupo == 1 ~ "Grupo 1",
    grupo == 2 ~ "Grupo 2",
    grupo == 3 ~ "Grupo 3",
    grupo == 4 ~ "Grupo 4",
    grupo == 5 ~ "Grupo 5",
    TRUE ~ "Sin muestra representativa")
)

colores <- mapa_grupos %>%
  dplyr::select(-nombre_barrio) %>% 
  distinct(color, .keep_all = TRUE) %>% 
  arrange(desc(grupo))


leaflet(data = political) %>% 
  addTiles() %>% 
  addProviderTiles(providers$OpenStreetMap) %>% 
  addPolygons(fill = TRUE, stroke = TRUE, weight = 2, color = mapa_grupos$color, 
              label = as.character(political$NOMBRE),
              popup = as.character(political$NOMBRE)) %>% 
  addLegend("bottomright", colors = colores$color, labels = as.character(colores$descripcion))
grupo_jerar_1<-grupos_barrios[grupos_barrios[,2]==1,]
grupo_jerar_2<-grupos_barrios[grupos_barrios[,2]==2,]
grupo_jerar_3<-grupos_barrios[grupos_barrios[,2]==3,]
grupo_jerar_4<-grupos_barrios[grupos_barrios[,2]==4,]
grupo_jerar_5<-grupos_barrios[grupos_barrios[,2]==5,]
nrow(grupo_jerar_1)
## [1] 271
nrow(grupo_jerar_2)
## [1] 26
nrow(grupo_jerar_3)
## [1] 9
nrow(grupo_jerar_4)
## [1] 1
nrow(grupo_jerar_5)
## [1] 1